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Abstract 

We are presenting a simple and numerical stable algorithm for the solution of the cone projection 
problem which is suitable for relative small data sets and for simulation purposes needed for con- 
vexity tests. Not even one pseudo-inverse matrix is computed because of a proper Gram-Schmidt 
orthonormalization process that is used. 
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THE CONE PROJECTION PROBLEM 



to observe that we have strict inequalities: 



Wc have the data set {xi, (j>i),i = 1,2, . . . ,n which 
has emerged from a convex function f at least C*-^-* [xi, Xn] 
by the process: 



Xi < X2 < ■ ■ ■ < Xn 



<Pi = f{xi) + Ei, iid{0, /„) 



(1) 



We want to find the vector y that has the smallest 
euclidcan distance from subject to the requirement of 
convexity Ay > 0, thus we have to solve the next primal 
optimization problem: 



so starting from the definition of convexity we proceed 
to the inequalities: 



Vi+2-Vi+l Vi+l-Vi 
a;i_|.2 — — Xi+i — Xi 



(2) 



subject to: {—Ay) < 



There are two equivalent versions for the matrix A of 
the convexity inequalities constraints. The first one is is 



iyi+2 - Vt+i) {xi+i - Xi) > {yi+i - yi) {xi+2 - Xi+i) 

{xi+2 - Xi+i) yi + {xi - Xi+2) Vi+i + {xi+i ~ Xi) yi+2 > 

(3) 

By constructing now all the above inequalities for 
i — 1,2, . . . ,n ~ 2 wc have formulated the matrix A'*^. 
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X3 - X2 



(4) 



Xn — 1 Xn — 2 Xn ^n— 1 -^n — 2/ 



The second way is obtained if we have equal spaced again with i = 1,2, . . . ,n — 2 and create the matrix A^'" 
Xi-data. Then it is easy to eliminate the same positive 



quantity Aa; ~ I'j+i — Xj from all inequalities: 

{xi+2 - Xi+i) yi + (xi - Xi+2) yi+i + {xi+i - Xi) yi+2 > 
{Ax) - {2Ax) y,+i + {Ax) y,+2 > 
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(6) 



yi - 2y^+i + y,+2 > 



Lemma .1. The polar component p* of the vector de- 
(5) composition (f> = y* + p* with y* the solution of problem 



2 



is a linear combination of the negative rows of matrix 
A, while the coefficients are either zero (if the correspond- 
ing constraint is not binding - inactive) or positive (if the 
relevant constraint is binding - active). 

Proof. The Lagrangian function of the problem and the 
first order condition for y can be written as: 



L (y, A) = (2/ 



(y ^ <j>) + i-Ay) (7) 



dLiy,\) 



2 (y - (A) + {-A^) A = 



2 (2/ - 0) + i-A'^) A = 
or 

y^^+^A^X 



(8) 



So, for the optimal solution {y* , A*} we have that it also 
holds: 



or 
or 

T\ A* 



0-y 



or 



-An X* 



Thus the representation of the polar component of the 
data vector, following the definitions of Q and i in the 
basis of the negative rows of A is half the Lagrange coef- 
ficient vector of the optimization problem [2] The coeffi- 
cients are zero or positive if the corresponding constraint 
is inactive or active respectively, due to Karush Kuhn 
Tucker complementarity slackness conditions. □ 



A NUMERICAL STABLE GEOMETRIC 
ALGORITHM FOR CONE PROJECTION 

An illustrative example 

Example .1. Let' s start with a common convex func- 
tion: 



f{x) = x\xe[0,i] 



Let the vectors a; ~ (0, ^, 1, |, 2) and < 



(10) 
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2 ' ' 2 

as presented at Figure [7] where we have drawn also the 
chord connecting {xi,(/)i) and (x^^ip^). If our data was 
convex then all {xi,(j)i) should lie above the chord, so 
clearly we have not convexity here. 

The demand for convexity takes the form of the next 
inequality constraints, using matrix A^"^ because of the 



f(x) = x'' (dot) 
chord ( solid line ) , errored data ( circle - cross ) 




FIG. 1: The statement of the convex projection problem for 
n = 5 
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equal 



A: 



Ay>0 
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(11) 



We define the matrix: 

-12-10 

i? = = I -1 2 -1 

0-12 -1, 

We pre-multiply vector cj) by R: 




If all components of the result were negative, then for the 
matrix A it should hold A(f> > 0, so our data should be 
convex. So, if we seek for the greatest deviation from 
convexity then it is natural to pick the component that is 
the greatest positive. This is compatible with (i) deviation 
from convexity and with (ii) Lemma\7T\ 
Here we observe the greatest entrance to be the 3rd one, 
so we pick up the 3rd row of R as the best chance to 
obtain a component of the polar vector and proceed by 
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taking the orthogonal projection of (f> onto that row: 
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Pi = Ri ^1 



yi= (j)- Pi 
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Now the greatest entry is the 2nd one, so we pick up the 
2nd row of R and continue by taking the matrix with the 
2nd and 3rd rows of R. It is important to notice that we 
are always sorting our indices in ascending order. 

/O 0\ 

-1 
2 -1 
-1 2 

yo -ij 

Now there exist two ways for projecting our data (f> on 
the two columns of R2 ■' 

1. The traditional way, i.e. the OLS estimator, which 
involves the pseudo inverse matrix and implies 
many numerical instabilities 

2. The new proposed way of taking the projection on 
the orthonormal base produced from them via the 
Gram-Schmidt procedure. 

We choose 2'"^ way and first construct for R2 with Gram- 
Schmidt the matrix with orthonormal columns: 
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Then we just take the projections of (p on the two 
columns of V: 



fi2 = V^<l> 









< 


\ 3730 1 




\ 20 / 





0.3061862179 
0.8215838362 



The reader who is familiar with the active set methodol- 
ogy has to notice that our vector is not identical anymore 
to the X vector of that method, because of the orthonor- 
malization process. This is the cost for the numerical 
stabilization of our algorithm. We continue executing our 
algorithm: 



P2 = V P2 



y2 




R2V2 = 



We observe that there exist no polar edge vector to be 
inserted in our algorithm, so we exit with the solutions: 



( \ 
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P = P2 
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We check again our solution: 
Ay* = 



It is expected to find that we will have two distinct lines 
for our cone projection plot and the second line has to be 
the OLS line for the set {{xi,(j)i),i = 2,3,4,5}, because 
only then we have a sequence of vanishing constraints. 
This fact is easily observed at Figure\^ 



The algorithm 

By increasing the dimension of our problem until 
a rather big n we continue to apply the same actions, 
i.e. we have estabhshed an algorithm for cone projection. 

We start by testing if our data is convex, so there is no 
need for cone projection at all. If it is not convex, then 
we multiply it by the R matrix and seek for the maximum 
component and for its position. That direction is more 
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f(x) = (dot) 
chord, errored data ( circle - cross ) 
cone projection (box) , polar vector ( diamond 




FIG. 2: The geometry of convex projection problem for n — 5 
in xy-plot 



probable to be an edge of the polar cone, so we find the 
projection of our data onto the i*'' row of R matrix. Now 
we have found the first approximations of the vectors p 
and y = (j) — p. We multiply again this y with R (b ~ Ry) 
and seek again for the maximum component and for its 
position. The new direction forms a set together with the 
previous one and we always sort the indices. The sorted 
indices construct the X matrix by taking the correspond- 
ing rows of R matrix as the columns of X. Then we apply 
the Gram-Schmidt orthonormalization procedure on the 
columns of X and construct the matrix V. Now our p 
vector can be calculated and then we find the next p and 
y approximation. We continue our algorithm until we 
reach at least one of the next three termination criteria: 

1. The algorithm is terminated if some b vector is 
'practically' zero 

2. The algorithm is terminated if there is no improve- 
ment in the value of b 

3. The algorithm is terminated if next index i of 
i?— row has already been chosen. 

Finally we exit from the algorithm with the set of in- 
dices J, where we have that the convexity constraints are 
satisfied as equalities (the active set indices, but without 
calculating the corresponding Lagrange coefficients), the 
polar component p* and the cone projection component 
y* of our initial data 0. We do not compute even one 
time any kind of pseudo-inverse matrix, which is the 
fundamental tool of every regression technique. This 
is due to the use of Gram-Schmidt orthonormalization 
process in order to do our orthogonal projections. 



This makes the algorithm numerical stable for using 
it for simulation purposes: we can establish the cone 
projection solution for every random set of vectors. 
This cannot be done with the traditional OLS solution, 
because of the existence of almost singular matrices for 
floating point arithmetic computations. The pseudo- 
code of the Algorithm is presented below. 

A Gram-Schmidt polar basis 
Cone Projection Algorithm 

Find y* ~ argmin \\y — (/)||2 



subject to Ay > 



INITIALIZE 

{61,62, J = {},i?, 



-A,b = Rcb, 6o, 



Id 



>0} 



• IF {6>0} THEN {p = 0,y^(f)} 
BREAK 

• ELSE 




FIRST PROJECTION 

= arg (bj = s) J = 

i=l,...,n-2 

p ~ projectnCp 
y = (j) - p,b = Ry 
NEXT PROJECTIONS 



sort (J U {i}) 



max bj i 



i=l,...,n-2 



arg {bj = s) J = sort (J U {i}) 

j=l,...,n-2 



IF {s < ei} THEN BREAK 

DO WHILE \\b-boid\\, > £2 



r 




V ~ GramSchmidt{X) J 

p ~ projectv4' 
y = (f> - p,b ^ Ry 

- arg {bj = s) J = sort (J U {i}) 

j=l,...,n-2 

- IF {s > ei} THEN 

{J = sort{JU{i})} 

-IF {ie J} THEN BREAK 
ELSE BREAK 

END DO 

CHECK SOLUTION \{y,p)\ < ei 
RETURN {J,p,y} 

We have developed our algorithm in four proper 
languages: 

1. First in Maple symbolic algebra system, where with 
just one page code we are able to execute our algo- 
rithm in absolute accuracy by using rational num- 
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bers as input data. 

2. Second in R suite, where we have floating point 
arithmetic, but it is an 'aher ego' for the statistician 
community. 

3. Third in Matlab/Octave, for those who are familiar 
with the benefits of them. 

4. Fourth in FORTRAN, one of the fastest ways to 
execute any numerical algorithm. 

CONCLUSION 

The presented algorithm for solving the cone projec- 
tion problem is quite simple because: 

• We don' t take care about the sign of Lagrange 
multipliers since we don' t compute them 

• We just include one component of the polar basis 
every time 

The algorithm is numerical stable for every kind of initial 
random vector (j) because all projections are done via a 



Gram-Schmidt procedure and not with the common OLS 
pseudo-inverse matrix, which is very often close to singu- 
lar for floating point arithmetic computations. 
The algorithm is useful for convexity tests where we need 
to compute the weights of the weighted or Beta dis- 
tribution that emerges for the corresponding statistical 
test, see for example [Ij. 
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